Is graphene in vacuum an insulator? 
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We present evidence, from Lattice Monte Carlo simulations of the phase diagram of graphene as 
a function of the Coulomb coupling between quasiparticles, that graphene in vacuum is likely to be 
an insulator. We find a semimetal-insulator transition at a"'* = 1.11 ± 0.06, where Og ~ 2.16 in 
vacuum, and Ug ~ 0.79 on a Si02 substrate. Our analysis uses the logarithmic derivative of the 
order parameter, supplemented by an equation of state. The insulating phase disappears above a 
critical number of four-component fermion flavors 4 < Nj'^^^ < 6. Our data are consistent with a 
second-order transition. 
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Graphene, a carbon allotropc with a two-dimensional 
honeycomb structure, has become an important player 
at the forefront of condensed matter physics, drawing the 
attention of theorists and experimentalists alike due to its 
challenging nature as a many-body problem, its unusual 
electronic properties and possible technological applica- 
tions (see Refs. [l|, [^ and references therein). Graphene 
also belongs to a large class of planar condensed-matter 
systems, which includes other graphite-related materials 
as well as high-T^ superconductors. 

A distinctive feature of graphene is that its band struc- 
ture contains two degenerate 'Dirac points', in the vicin- 
ity of which the dispersion is linear, as in relativistic the- 
ories 3] . The low-energy excitations in graphene are thus 
Dirac quasiparticles of Fermi velocity v ~ c/300, where c 
is the speed of light in vacuum. These are described by 
the Euclidean action 
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(1) 



where g^ — c^/eq for graphene in vacuum, ^p^ is a 
four-component Dirac field in 2+1 dimensions, Aq is 



a Coulomb field in 3- 
graphene, and 



1 dimensions, Nf = 2 for real 



D[Ao] = 70(^0 + ^A,) + vjA. * = 1, 2 (2) 



where the Dirac matrices 7„ satisfy the Euclidean Clif 
ford algebra {7„,7^} = 26 



f,, luj -t,u- The strength of the 

Coulomb interaction is controlled (as can be shown by 
rescaling t and Aq) by a„ = e^/(47rueQ), which is the 
graphene analogue of the fine-structure constant a ~ 
1/137 of quantum electrodynamics (QED). 

Despite the similarities with QED, the smallness oiv/c 
in graphene has non-trivial consequences: Coulomb inter- 
actions between the quasiparticles are essentially instan- 
taneous, thus breaking relativistic invariance which is re- 
flected in Eq. ([1]). The analogue of the fine-structure con- 
stant a„ "^ 300 a in graphene, such that the low-energy 



properties resemble QED in a very strongly coupled 
regime. This provides an exciting opportunity for the 
study of strongly coupled theories, within a condensed- 
matter analogue that can be experimentally realized with 
modest equipment. 

Notably, Eq. ([1]) satisfies a chiral U(2A^^) symmetry 
which can break spontaneously at large enough Coulomb 
coupling, generating a gap in the quasiparticle spec- 
trum. Whether such an effect occurs in real graphene 
is an open issue from the experimental point of view (see 
however Ref. [J], where a substrate-induced gap is re- 
ported). On the theoretical side, dynamical gap genera- 
tion is described by a quantum phase transition due to 
the formation of particle-hole bound states. However, in 
such a strongly coupled regime, even a qualitative anal- 
ysis should be non-perturbative. This is especially true 
for graphene in vacuum, where a„ attains its maximum 
value, while it is partially screened in the presence of a 
substrate. While the semimetallic properties of graphene 
on a substrate are well established, the issue of whether 
a transition to an insulating phase occurs in the absence 
of a substrate remains unsettled. 

This problem has been studied using perturbative as 
well as non-perturbative methods [5|j la LDi ISl • The lat- 
ter, which are typically based on a gap equation, yield 
an infinite-order transition to an insulating phase above 
a critical coupling. On the other hand, the results of 
laige-N r analyses [5|, |9|, llO| find that Coulomb inter- 
actions flow towards a non-interacting flxed point un- 
der renormalization-group (RG) transformations, being 
therefore unable to induce a transition. However, a num- 
ber of uncontrolled approximations are involved, such as 
the reliance on large-A'^^ results that may break down 
for small N ^ , various approximate treatments of the gap 
equation kernel, as well as the linearization of the result- 
ing integral equation (see Appendix B of Ref. [8|). 

We set out to characterize the phase diagram of 
graphene in a controlled fashion, which entails a lattice 
Monte Carlo approach and analysis of the chiral conden- 
sate, which is the order parameter for a transition into an 



insulating charge-density-wave phase. Such an approach 
is non-perturbative, takes full account of quantum fluc- 
tuations, and has been extensively used [ll|, [ij, llj, ll^l 
in the study of 2-hl (QEDg) and 3-^1 (QED4) dimen- 
sionalQED, but not for studies of graphene (see however 
Ref. [l5|, where a model for the strong-coupling limit is 
investigated). 

To this end, we discretize the pure gauge part of Eq. ^ 
according to (for a recent overview, see Ref. [l6() 
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where the (dimensionless) lattice coupling f3 = f/5^, Oq 
is the lattice gauge potential, n = (rig, . . . , rig) denotes a 
site on the space-time lattice, and e is a unit vector in 
the direction fi. For studies of chiral phase transitions, 
staggered fermions [17| are a preferred choice, since chi- 
ral symmetry is then partially preserved. As N stag- 
gered flavors correspond to Nf = 2N continuum Dirac 
flavors [18[, it suffices (for real graphene) to set iV = 1, 
which gives 
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FIG. 1: (Color online) Chiral condensate a for Nj —2,4,6 as 
a function of /3 and mo, with lines intended to guide the eye. 
The lattices are of extent L"^ x L^, such that the fermions live 
in a 2-1-1 dimensional cube of size L, while the gauge bosons 
also propagate in the z-direction of length L^. For each /3, 
the value of L is given next to the datapoints. All results are 
for L^ = 8, as larger values had no discernible effects. For 
each datapoint ~ 300 uncorrelated gauge configurations were 
generated. The statistical uncertainties, which are compara- 
ble to the size of the symbols, were obtained by the jackknife 
method [24|. Finite volume effects are largest for small /3. 



where the Xn ^^^ staggered fermion spinors, and (m, n) 
are restricted to a 2-1-1 dimensional sublattice. The 
invariance of Eq. ([1]) under spatially uniform, time- 
dependent gauge transformations is retained by coupling 
the fermions to the gauge field via C/q = exp(i6'Q). The 
staggered form of D is 
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where 77^ „ = (— l)"o and r]2^n = (— l)"o+"i. The mass 
term breaks chiral symmetry explicitly, generating a non- 
zero condensate which is otherwise not possible at finite 
volume. Extrapolation to ttt-q = is thus required. Upon 
integration of the fermionic degrees of freedom, the path 
integral is governed by the effective action 



SMOo] - -NlndetiD^^jUo]) + SUdol 



(6) 



such that P[Oq] = exp{~ScB[OQ]) defines the Monte Carlo 
probability measure, ft is straightforward to show that 
the determinant is positive definite. We have sampled 
P[9q] using the Metropolis algorithm, updating 9q at ran- 
dom locations and evaluating the fermion determinant 
exactly. Our approach has been tested against known 
results for QEDg and QED4. The data of Ref. [lH on 
the chiral condensate of QEDg have been accurately re- 
produced for multiple values of Nr, along with several 
randomly chosen datapoints from Refs. []J 
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The chiral \J(2Ni:) symmetry of the continuum theory 
can only partially be realized on the lattice if the doubling 
problem is to be avoided [19]. In particular, only a global 
\J{N)x\J{N) symmetry remains upon discretization 18|. 
We focus on the spontaneous breakdown of this symme- 
try to a V{N) subgroup, characterized by a condensate 
a = (xx) 7^ in the limit TOq — > 0, which marks the 
appearance of a gap in the quasiparticle spectrum. Our 
results for a are presented in Fig. [T] for /3 = 0.05, . . . , 0.5 
and TTiQ = 0.010, . . . , 0.025 (in lattice units). 

Our data for Nr = 2 in Fig. [1] are suggestive of a crit- 
ical coupling /3c ~ 0.06 . . .0.09, below which a survives 
in the limit ttiq — > 0. More significantly, the suscepti- 
bility Xi = da/drriQ shown in the right panel of Fig. [2] 
exhibits a maximum which tends towards P^ as rrig is de- 
creased. The (much more limited) data for Nf = 4 (not 
shown) have a similar maximum around /3 ~ 0.03. As N^ 
is increased, a obviously becomes suppressed, vanishing 
between N^ — 4 and Nf — 6. This agrees with large- A^^: 
results jlOJ I that yield a quantum critical point (and there- 
fore no condensate) in the limit /? — ^ 0, and is consistent 
with recent Monte Carlo studies of that limit [15]. Our 
results for small A^^ establish that the strong-coupling 



critical point disappears below A^S''* , 



with 4 < iV^"* < 6. 



For a quantitative determination of (3^, we compute 
the logarithmic derivative R oi a with respect to tuq, 



R = 



dim 



d In TOp 



da 
drrir, 



(7) 
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FIG. 2: (Color online) Left panel: Logarithmic derivative 
i? as a function of mo for different /3. Right panel: Chiral 
susceptibility Xi as a function of (3 for different rrio. All data 
are for Nf = 2, with lattice sizes identical to those of Fig. [l] 
The lines are intended as a guide to the eye. 



according to the method of Ref. [20] . 
the behavior of R is as follows: R -^ 



In the limit : 



0, 



1 in the chirally sym- 



metric (unbroken) phase, where a oc TOq. At the critical 
couphng j3 — fi^ one finds that R -^ 1/5, where (5 is a 
universal critical exponent. R vanishes in the sponta- 
neously broken phase, where cr 7^ for ttiq -^ 0. The 
data in Fig. [5] (left panel) indicate that chiral symme- 
try is spontaneously broken for (3 — 0.067, but remains 
unbroken for (5 — 0.077, from which we conclude that 
/3^ = 0.072 ± 0.005. This estimate can be refined by use 
of larger lattice volumes and smaller values of TOq . 

A more precise determination of [3^ requires an equa- 
tion of state (EOS) of the form m^ = f{a,f3) for the 
extrapolation TOq — > 0. We have considered the EOS 
successfully applied [1^, ll3[ to QED4 



^4' 



m,X{f]) = y(/3)/i(a) + /3(a). 



(8) 



where X{l3) and Y{P) are expanded around f3^ such that 
X(/3) = Xo+X,{l-P/f3,) andy(/3) = Y,{l-P/f3,). The 
dependence on a is given by /^(cr) = a'' and fr^{cr) — ct*, 
which allows for non-classical critical exponents S and 
^ [11 m, where b=S- 1/p. A x^ fit of Eq. ^ to the 



data in Fig. [T] is given in Fig. [3] (left panel), along with a 
"Fisher plot" of cr^ versus ra^/a (right panel). Deviations 
from classical mean-field behavior manifest themselves in 
the Fisher plot as curvature in the lines of constant (3. 
The mostly straight lines in Fig. [3] (right panel) indicate 
that the data should be well described, up to finite-size 
effects, by the classical values ^ = 3 and 6=1, which is 
confirmed in Fig. [3] (left panel) . 

For all values of ttiq, the finite- volume effects appear 
to be dominated by a dynamically generated correlation 
length in the region that we identify as the spontaneously 
broken phase. To gauge the impact of such effects on /3^, 



several fits were performed from L = 8 up to the main 
fit with L = 16. Our analysis indicates that /3^ is sta- 
ble around ~ 0.075, even though the datapoints at the 
smallest /3 shift due to finite-volume effects. In contrast, 
for N f = 6 chiral symmetry remains unbroken, and thus 
finite- volume effects remain small in the limit /3 — > 0. 
The stability of j3^ can be understood in terms of the 
results for R in Fig. [21 since any fit to the condensate 
should be consistent with the susceptibility as well, and 
R combines both pieces of information. The character of 
R clearly changes around (3 = 0.075, in direct correlation 
with the fitted values of (3^. If one allows for extreme 
modifications, such as non-classical critical exponents, or 
forces the fit to account for all the datapoints, values of 
(3^ as low as 0.060 may be found, although at the price of 
a much worse fit to the data. By inclusion and exclusion 
of different sets of datapoints, a realistic (though some- 
what model-dependent) estimate of the critical coupling 
is (3^ — 0.0755±0.0008. Possible deviations from classical 
mean-field behavior with fi=cr and f^ = a^ are below 
the resolution of the present study. 

The results in Fig. [3] suggest a second-order transi- 
tion with classical exponents, unlike Refs. [a, Isl], where 
an infinite-order transition was found. In this situation, 
further investigation is clearly called for. While the sen- 
sitivity of our analysis increases for smaller ttiq, larger 
lattice volumes are also required to keep finite- volume ef- 
fects under control. A similar EOS analysis has recently 
been performed by Hands and Strouthos (Ref. 15]) for a 
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FIG. 3: (Color online) Left panel: x^ fit to the data of Fig. [T] 
and extrapolation to rrio — for A^^ = 2 using Eq. ((8} with 
X(),Xi,Yi and (3^ as free parameters. The points with largest 
finite- volume effects have been excluded from the fit. The 
optimal parameter values are /3c — 0.0755 ± 0.0003, Xq = 
0.195 ±0.003, Xl = -0.089 ±0.001 and Yi = -0.091 ±0.001. 
The uncertainties are purely statistical. Right panel: Fisher 
plot of a^ versus rrio/cr for the data of Fig. [T]with Nf = 2. The 
lines connect datapoints with identical /3, such that straight 
lines indicate mean-field behavior according to Eq. ((8)1 . At 
/3c, the extrapolation crosses the origin. 



graphene-like theory with a zero-range interaction. Un- 
fortunately, a meaningful comparison is not possible at 
this time, as their parameter \l g^ cannot be identified 
with our gauge coupling /3, except in the strong-coupling 
limit /3 ^ 0. 

A comparison with experiment necessitates a discus- 
sion of renormalized quantities. While test charges re- 
main unscreened as the fluctuations of the fermion ac- 
tion are confined to 2 + 1 dimensions, the physical value 
of fi^ may be affected by renormalization of v due to the 
breakingof relativistic invariance. Large- TV j- results sug- 
gest [5|, |lO[ that the Coulomb interaction renormalizes v 
logarithmically toward larger values, thereby decreasing 
a„ slightly from the bare value, and strengthening our 
conclusions for graphene in vacuum. Available experi- 
mental evidence [J, |6| indicates that velocity renormal- 
ization is at most a ^ 20% effect, and of phononic rather 
than Coulombic origin. 

Summarizing, we have found that graphene should be- 
come insulating at a critical coupling a^'''* = l/(47r/?j,) — 
1.11 ± 0.06, where /?c = 0.072 ± 0.005. This should be 
compared with a ~ 2.16 in vacuum, and a ~ 0.79 
on an Si02 substrate (using the experimental value v ~ 
10^ m/s). These findings are in line with the observed 
semimetallic properties [1] of graphene on a Si02 sub- 
strate, and predict that the Coulomb interaction in sus- 
pended graphene should induce a gap in the quasiparticle 
spectrum. Within the accuracy of the present study, the 
transition appears to be of second order. 

Ultimately, the observation of the insulating phase is 
dependent on the size of the induced band gap. How- 
ever, the prediction of a dimensionful observable requires 
the matching of a lattice quantity (other than the gap 
itself) to the corresponding experimental value. This ap- 
plies to finite temperature studies as well, where it is 
necessary to fix the absolute temperature scale. An in- 
triguing possibility is that the observed nanoscale ripples 
in suspended graphene ^21.] may provide the necessary 
information, as such corrugations can by described |22l | 
by means of external gauge fields with known dimension- 
ful properties. Exploratory work in this direction, along 
with a more accurate study of the transition properties, 
is in progress [2^ . 
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